!======================================================================================
integer function itest_mesh()

! 1) Testing for remeshing
! find the smallest angle in each 4 subtriangles in each qudralateral   
! and compare to angle_for_remeshing
 
! 2) Find the zone where it is needed to be remeshed:
! go from the left to the right and backwards
!
! G.Ito if ny_inject=3, then only remesh during tectonic periods
!=====================================================================================
include 'precision.inc'
include 'params.inc'
include 'arrays.inc'
dimension iv(4),jv(4),angle(3)


itest_mesh = 0
iac_rem=0				!accretion check deactivated G. Ito

! if remeshing with adding material on the sides then
! remesh at pre-given shortening
! dx_rem*dx - critical distance of shortnening
if( mode_rem .eq. 11 ) then
    testcr = dx_rem * rxbo / (nx-1)
    shortening = abs(cord(1,nx,1) - cord(1,1,1) - rxbo)
    if ( shortening .gt. testcr ) then
        if( dtout_screen .ne. 0 ) then
            print *, 'Remeshing due to shortening required: ', shortening
        else
            call SysMsg('TEST_MESH: Remeshing due to shortening required')
        endif
        itest_mesh = 1
        return
    endif
end if

!test by dx in the first layer
dxratmax = 5.
if(ny_inject.ge.1) dxratmax = 8.		!G.Ito was if (ny_inject.eq.2)
dxmin = 1.d20
dxmax = -1.d20
do i = 1, nx-1
    dx = cord(1,i+1,1) - cord(1,i,1)
    if( dx .lt. dxmin ) dxmin = dx
    if( dx .gt. dxmax ) dxmax = dx
enddo

if( dxmin.lt.0 .or. dxmax/dxmin.gt.dxratmax ) then
!  if ((ny_inject.ne.3) .or. (rateM_mech.ge.1.0)) then			!G.Ito 12/06
    if( dtout_screen .ne. 0 ) then
        write(*,'(a35,2f10.3)'), 'Remeshing due dxmax/dxmin criteria:', dxmin/1.e3, dxmax/1.e3
    else
        call SysMsg('TEST_MESH: Remeshing due dxmax/dxmin criteria required')
    endif
    itest_mesh = 1
    return
!  elseif (tinj.lt.(ptect-50.0*sec_year)) then  
!    if( dtout_screen .ne. 0 ) then
!        write(*,'(a35,2f10.3)'), 'Remeshing due dxmax/dxmin criteria:', dxmin/1.e3, dxmax/1.e3
!    else
!        call SysMsg('TEST_MESH: Remeshing due dxmax/dxmin criteria required')
!    endif
!    itest_mesh = 1
!    return
!  endif  
endif

! This check deactivated by G. Ito.  (it wasn't working correctly
! for variable grid spacing because dx_int=abs(cord(1,2,1)-cord(1,1,1))
! test mesh in the case of accretion at the center
!if (ny_inject.eq.2) then
!check the width of the two elements at the center
!dx_accr = abs(cord(1,iinj-1,1)-cord(1,iinj,1))
!if (dx_accr.ge.1.05*abs(dx_init)) then
!if (dx_accr.ge.2.*abs(dx_init)) then
!itest_mesh = 1
!iac_rem = 1 
!endif
!    if( dtout_screen .ne. 0.and.itest_mesh.eq.1 ) then
!        print *, 'Remeshing due to accretion required.'
!        mloop = 0
!    endif
!endif
! test by the angle
if( mod(nloop,ntest_rem).ne.0 ) return

pi = 3.14159265358979323846
degrad = pi/180.
raddeg = 180./pi
anglemint = 180. 
imint = 0
jmint = 0

do 1 i = 1, nx-1
    do 1 j = 1,nz-1

        ! loop for each 4 sub-triangles
        do ii = 1,4
            if (ii.eq.1) then
                iv(1) = i ; jv(1) = j ; iv(2) = i ; jv(2) = j+1 ; iv(3) = i+1 ; jv(3) = j
            elseif (ii.eq.2) then
                iv(1) = i ; jv(1) = j+1 ; iv(2) = i+1 ; jv(2) = j+1 ; iv(3) = i+1 ; jv(3) = j
            elseif (ii.eq.3) then
                iv(1) = i ; jv(1) = j ; iv(2) = i ; jv(2) = j+1 ; iv(3) = i+1 ; jv(3) = j+1
            elseif (ii.eq.4) then
                iv(1) = i ; jv(1) = j ; iv(2) = i+1 ; jv(2) = j+1 ; iv(3) = i+1 ; jv(3) = j
            endif
            iv(4) = iv(1)
            jv(4) = jv(1)

            ! Find all angles using vector dot product a*b = |a||b|cos(a)
            do k = 2,3
                xa = cord(jv(k+1),iv(k+1),1)-cord(jv(k),iv(k),1)
                ya = cord(jv(k+1),iv(k+1),2)-cord(jv(k),iv(k),2)
                xxal = sqrt(xa*xa + ya*ya) 
                xb = cord(jv(k-1),iv(k-1),1)-cord(jv(k),iv(k),1)
                yb = cord(jv(k-1),iv(k-1),2)-cord(jv(k),iv(k),2)
                xxbl = sqrt(xb*xb + yb*yb) 

                angle(k)=raddeg*acos((xa*xb+ya*yb)/(xxal*xxbl))
            end do
            angle (1) = 180.-angle(2)-angle(3)

            ! min angle in one trianle
            anglemin1=min(angle(1),angle(2),angle(3))  

            ! min angle in the whole mesh
            if( anglemin1 .lt. anglemint ) then
                anglemint = anglemin1
                imint = i
                jmint = j
            endif

        end do

1   continue
2   continue

if( dtout_screen .ne. 0 ) &
    write (6,'(A,F5.2,A,I3,A,I3)') '        min.angle=',anglemint,' j=', jmint, ' i=',imint

! check if the angle is smaller than angle of remeshing  
! Check for ny_inject=3 & force regridding during high-stress state of tectonic cycle.  G. Ito 12/11/06
if (anglemint .le. angle_rem) then
!     if ((ny_inject.ne.3) .or. (rateM_mech.ge.1.0)) then
      if( dtout_screen .ne. 0 ) then
          write(6,'(a54,2f12.4)'), 'Remeshing due to angle required.  time (kyr), angmin =', &
	  time/sec_year/1.e+3, anglemint
          mloop = 0
      else
          call SysMsg('TEST_MESH: Remeshing due to angle required.')
      endif
      itest_mesh = 1
      iac_rem = 0
!     elseif (tinj.lt.(ptect-50.0*sec_year)) then  
    			!Need ~50 yrs to re-equilibrate after regridding
!       if( dtout_screen .ne. 0 ) then
!           write(6,'(a54,2f12.4)'), 'Remeshing due to angle required.  time (kyr), angmin =', &
! 	  time/sec_year/1.e+3, anglemint
!           mloop = 0
!       else
!           call SysMsg('TEST_MESH: Remeshing due to angle required.')
!       endif
!       itest_mesh = 1
!       iac_rem = 0
!     endif
endif

return
end
